program PS

        ! "use" is the inheritance mechanism for F90:
        use mConstants
        use mObject
        use mUniverse
        use mLPS
        use mPSFittingFunctions
        use mMisc

        implicit none

        real(D) :: smithp3d, linearp3d, pdp3d
        real(D) :: k(200), kfactor, l1(200), smith1(200), pd1(200)
        real(D) :: sig8,dsig8,kr,Wsig,dlnk,twopisq
        real(D) :: om_m,om_v,om_b,h,A_s,n,w,wa,theta,ll
        integer :: i,loop,nloop

        nloop = 100
        twopisq = 19.7392_D

        kfactor = exp(log(1.e10_D)/199._D)
        k(1) = 1.e-5_D
        do i = 2, 200
          k(i) = k(i-1) * kfactor
        enddo
        
        open(10,file='WWaPrior.dat',status='old')
        do loop = 1,nloop
          print *,'loop = ',loop
          read(10,*) om_b,om_m,theta,n,wa,w,A_s,h,ll
          om_b = om_b / h / h
          om_m = om_m / h / h
          om_m = om_m + om_b
          om_v = 1.0_D - om_m

          call def_universe(om_m,om_v,om_b,h,A_s,n,w,wa)

          do i = 1, 200
            smith1(i) = smithp3d(k(i),1._D);
            l1(i) = linearp3d(k(i),1._D);
            l1(i) = l1(i) * k(i)**3 /twopisq
            pd1(i) = pdp3d(k(i),1._D);
          enddo

          if (loop .eq. nloop) then

            write( 6,*) '     k        Lin       Smith        PD'

            sig8 = 0._D
            do i = 1, 200
              kr = k(i)*(8.0_D)/h
              Wsig = 3._D*(sin(kr)-kr*cos(kr))
              dlnk = log(kfactor)
    
              dsig8 = W*W*l1(i)*dlnk
              sig8 = sig8 + dsig8
              write( 6, fmt = '(5e11.3)' ) k(i),l1(i),smith1(i),pd1(i), pd1(i)/smith1(i)
            enddo
            write(6,*) 'sigma_8 = ',sig8
          endif
        enddo

end program PS
